Comparison of Perron and Floquet eigenvalues in age 
structured cell division cycle models 



"oo ■ Jean Clairambault"'^, Stephane Gaubert^'^ and Thomas Lepoutre^ 

o ■ 
o 

(N 

O ■ " INRIA, projet BANG, Domaine de Voluceau, BP 105, 78156 Le Chesnay Cedex France 

^ '. ^ INSERM U 776, Hopital Paul-Brousse, 14, Av. Paul-Vaillant-Couturier F94807 Villejuif cedex 

^ : " INRIA Saclay - Ile-de-France, projet MAXPLUS 

^ ' CMAP, Ecole Polytechnique, 91 128 Palaiseau Cedex, France 

_ ^ UPMC Univ Paris 06, UMR 7598, Laboratoire Jacques-Louis Lions, F-75005, Paris, France 

CM ■ 

< 

<^ • Abstract 

We study the growth rate of a cell population that follows an age-structured PDE with 
time-periodic coefficients. Our motivation comes from the comparison between experimental 
^ ■ tumor growth curves in mice endowed with intact or disrupted circadian clocks, known to exert 

i their influence on the cell division cycle. We compare the growth rate of the model controlled 

by a time-periodic control on its coefficients with the growth rate of stationary models of the 
Q ■ same nature, but with averaged coefficients. We firstly derive a delay differential equation 

^ . which allows us to prove several inequalities and equalities on the growth rates. We also 

discuss about the necessity to take into account the structure of the cell division cycle for 
chronotherapy modeling. Numerical simulations illustrate the results. 
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1 Cell cycle control and circadian rhythms 

The cell division cycle is the process by which the eukaryotic cell duplicates its DNA content 
and then divides itself in two daughter cells. This process is normally controlled by vaiious 
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physiological mechanisms that ensure homeostasis of healthy tissues, that control genome 
integrity (e.g. cyclins and cdks, p53, repair enzymes, etc.), launching programmed cell death 
(apoptosis) if the DNA is irreversibly damaged (see ||2TI for a complete presentation). The 
system of control has been extensively studied and modeled (see e.g. |[T4l [T6l l23l or |[26l ) 
using ordinary differential equations. The cell division can be modeled through branching 
processes (see 111), integral equations, delay differential equations (see H) and also many 
structured PDE models (for an overview, see HI [21 [191) where the structuring variables can be 
age ( |[22l ). size ( |[24l ) or more recently cyclin content (|[5ll6l[Tn). 

Most living organisms exhibit circadian rhythms (from Latin circa diem, "roughly a day") 
which allow them to adapt to an environment that varies with a periodicity of 24h. These 
rhythms can be observed even in the smallest biological functional unit, the cell. The problem 
we are studying is the growth of cell populations (undergoing the cell division cycle described 
above) under the pressure of circadian rhythms. Circadian rhythm effects on the cell cycle turn 
out to be important in tumor proliferation. This is observed by several experiments involving 
a major disruption of circadian rhythms in mice. In these experiments it can be seen that the 
growth of tumors is significantly enhanced in mice in which the pacemaker circadian clock 
has been drastically perturbed, either through neurosurgery, or through light-dark cycle dis- 
ruption (see e.g. llT3l[T2l '). Moreover, in the clinic, taking advantage of the influence exerted 
by circadian clocks on anticancer drug metabolism and on the cell division cycle has led in the 
past 15 years to successful applications in the chronotherapy of cancers, particularly colorectal 
cancer (see [;18] ). This motivates modeling the circadian rhythm in simple cell cycle models 
and studying these effects on the growth rate of a cell population. 




Figure 1 : Effects of the perturbation of light-dark cycle on tumor proliferation (reproduced from 
|[T2l ). In clock-perturbed mice (black dots), the tumor proliferates much faster than in control mice 
(white dots). (By courtesy of Elizabeth Filipski). 

Contrary to our first idea, the growth rate of a cell population described by a physiologically 
structured PDE model with time-periodic control is not necessarily lower than in a model of 
the same nature, but with a time-averaged control |[7l[8l[9l. 
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The goal here is twofold. Firstly we analyze how modeling assumptions lead to define 
various growth rates under the effects of circadian rhythms. Secondly we model the effect of 
chrono therapy on these growth rates. 

In the second section we recall the definition of these various growth rates, in terms of 
PeiTon and Floquet eigenvalues of a linear- Von Foerster- Mc-Kendrick model. We also discuss 
known inequalities between them. In the third section we study a simple division model, 
for which we establish (in Theorem |2ll strict inequalities comparing the growth rate in the 
stationary (Perron) and periodic (Floquet) cases. These inequalities are proved by studying a 
related time delay system (which is similar- to the one considered in t4j). This model is used 
to confirm the impossibility to derive a general comparison between the Perron and Floquet 
eigenvalues defined in the second section. In the fourth section, we give an argument for using 
multiphase models to represent chronotherapy, taking better into account the structure of the 
cell cycle and particularly the existence of various phases. We provide numerical simulations 
to illustrate our results. In a first appendix, we give the detailed proof of the existence of the 
solution of the eigenproblem, by applying the Krein-Rutman theorem. In a second appendix, 
we derive analytical formulse for the eigenelements in a specific multiphase case, which yield 
further information on their behavior and can be used to validate numerical experiments. 



2 The model 

2.1 The renewal equation 

We base our study on a cell population that foUows the classical renewal equation structured 
in age with periodic coefficients representing the effect of circadian rhythms 

^n(t, x) + -^n{t, x) + d{t, x)n{t, x) = 0, 
n{t, X = 0) = B{t, x)n{t, x)dx. 

Here n{t,x) represents the density of cells of age x in the cycle at time t, d{t,x),B{t,x) 
represent respectively the death rate, and the birth rate. Both these coefficients are T-periodic 
in time. We define the growth rate of the population in terms of an eigenproblem. The growth 
rate Xp {F for Floquet as for ODEs with periodic coefficients) is defined as the unique real 
number Xp, such that there is a solution to the problem 

' f 7V(t, x) + ^N{t, x) + [Xp + d{t, x)]N{t, x) = 0, 
< N{t, x = 0) = B{t, x)N{t, x)dx, (2) 
N > 0, T-periodic. 

We refer to ll20l for conditions of existence for A p (and to the appendix for the case of division 
models). 
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2.2 Comparison of eigenvalues 

We use the following notations. For a T-periodic function / we define, 
1 

(/) = — / f{t)dt the arithmetical average, 
Jo 

{f)g = exp J log f{t)dt^ the geometrical average, when / > 0. 

It may seem natural to introduce the following stationary problem (Perron eigenproblem), 
in which the death and birth rates are averaged 

f ^iVp(x) + [Ap + {d{x))]Npix) = 0, 

^P(O) = I^{B{x))Np{x)dx = 1, (3) 

_ Np{x) > 0. 

It is shown in HI 13 that, when B does not depend on time, the inequality Ap > Ap holds. 
In the present paper, we show that this inequality does not carry over to the case of a time 
dependent B. It should be noted, however, that there is a general inequality, established in 
Q, which relates Ap with the solution of the following eigenproblem in which an arithmetical 
average of the death rate is taken, whereas the geometrical average of the birth rate is taken, 

^N,ix) + [Xg + {d{x))]Ngix) = 0, 

^5(0) = I^{B{x))gNg{x)dx = 1, (4) 

Ng{x) > 0. 

Theorem 1 (171). The eigenvalues defined in ((21) and (I?]) satisfy 

Ap > Xg. 

This result suggests that there is no general inequahty between Ap and Ap, because the 
inequality which follows from convexity is Ap > Xg. Moreover, it follows from the standard 
arithmetico-geometrical inequality, 

Ap > Xg. 

Such a general comparison cannot hold between Ap and Ap, as shown in the next section. To 
go further we use a more specific model. 

3 A simple one-phase division model 
3.1 Model and main results 

We model the cell cycle with the following PDE which is a particular case of ([U, 

J |n(t, x) + ^n(t, x) + [d{t) + Ko'^{t)x[a,+oo[ix)Ht, x) = 0, 
\ n(t,0) = 2Ko7p{t) f^n{t,x)dx, 
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where Kq > is a constant, > is a T-periodic function with 

(V) = 1. (1) 

The term i^^oV'(*)X[a,+oo[ represents the division rate, d{t) is the apoptosis rate (we assume it to 
be T-periodic). We have denoted by xe the indicator function of set E. Finally, ip{t) represents 
a nonnegative periodic control exerted on division. As before we look for the growth rate A p 
of such a system. It is defined so that there is a solution to the Floquet eigenproblem, 

|7V(t, x) + ^N{t, x) + [Xf + d{t) + i^oV'(i)X[a,+oo[(^)] N{t, x) = 0, 
N{t, 0) = 2Ko^(t) XT ^(*' ^)^^' (2) 
iV > 0, T-periodic, 

and we normahze by 

rT /"OO 

/ / N{t,x)dxdt = 1. 
Jo Jo 

As we already know a general comparison result for the geometrical eigenvalue defined 
in (01), we are now only interested in the comparison of A^ and Ap, the latter quantity defined 
by requiring the existence of a solution to the PeiTon eigenproblem already defined in ([3]) which 
here reads 

' ^^P(^) + + (d) + KoX[a,+oo[{x)]Np{x) = 0, 
Np{0) = 2Kof;;^Np{x)dx, (3) 

_ iVp > 0, 
and we normahze Np by 



POO 

Np{0) = 2Ko / Np{x)dx = 1. 

J a 



We are interested in evaluating the effect of the periodic control ip{t) on the growth of the 
system. Therefore we denote by Ap(a,'i/') and by Ap(a) the above defined eigenelements so 
as to keep track of the problem parameters. 

The following theorem implies that there is no possible general comparison between Ap 
and Ap. 

Theorem 2. For all continuous positive T-periodic functions tp satisfying dU, we have 

Ap(a = T,^) =Ap(T) = Ap(a = T,l), (4) 
and for a in a neighborhood of T, we have, provided ip ^ 1 

Ap(a, ip) > Ap(a) = Ap(a, 1) for a < T, 

Ap(a, ip) < Ap(a) = Ap(a, 1) for a > T. 
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The proof of this theorem is presented in the next sections. The computations done in 
section 4. 1 insure that, without loss of generality, we can suppose d = 0. 

Numerical results are presented in figures |2] and |3] which illustrate this theorem. Graphi- 
cally, for fixed V'. this predicts firstly that the curves of Xpiajt/j) (Floquet curve) and Xp{a) 
(Perron curve) must cross each other for a = T, secondly that the Floquet curve should be 
above the Perron curve before (i.e., for a < T) the crossing and below this curve after it (i.e., 
for a > T). A possible interpretation is that for a better adaptation (in the sense of higher 
proliferation), the cell cycle should be shorter than 24h; an effect already observed in H. 

3.2 Proof of Theorem 2, part 1 (a delay differential equation) 

Throughout the proof, we use the shorter notations \p and Ap instead of Xpia, ip) and Xp{a) 
when there is no possible confusion. 

To find more information on Ai? we derive a delay differential equation. 
We integrate Q with respect to age over [a, cx3[. We get 

— / N{t, x)dx + N{t, oo) - N{t, a) + [Xp + Koip{t)] / N{t, x)dx = 0. 

"-t J a J a 

From the formula of characteristics and the boundary condition in 

N{t,a) = N{t-a,0)e-^P'', 

N{t, a) = 2Koe-^P''tp{t - a) N{t - a, x)dx. 

We set P{t) = N{t, x)dx. Since we have N{t^ oo) = (see the appendix) we obtain 
the delay differential equation 

P{t) + (Xp + i^oV'(i)) P{t) = 2Koijit - a)P{t - a)e-^^^ (5) 



3.3 Proof of Theorem 2, part 2 (equality of growth rates for a = T) 

The comparison between Ap and Ap is based on the following formula for Ap. 
Lemma 3. The Perron eigenvalue defined in ((21) satisfies 

Va > 0, ^^e^- = 1. (6) 

Proof. From ©, we have, for x > a, Np{x) = e-'^^p+^o)x+Koa_ -y^g ^^^^^ 
boundary condition and obtain 



/•oo 
J a 



1 = 2Ko e-^P". 

Xp + Ko 



□ 
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Corollary 4. The Perron eigenvalue defined in Q satisfies 

Va > 0, Ap > 0. 
Proof. This follows from Lemma|3]and the remark 

Va>0,VA<0, ^e^«<i. 



□ 



To obtain we divide ([5]) by P and find 

When we take the average over a period, we get (since P is T-periodic in time by its definition 
as is) 



= -{Xf + Ko) + IKoe-^^^l ^{t - a) 



Pit) 



Now we consider the particular case a = T. As P is T-periodic P{t — a) = P{t). Hence, for 
a = T, we arrive at 

^- + '^V-=/.*«-a)^\ = W) = l. (8) 



2K„ ' P(() 

This equality is the same for Ai? as the one described in lemma 1 for Ap. As we know that the 
mapping 

is increasing on [—Kq, +(X)[ from to +oo and is negative elsewhere, there is only one solution 
to ([8]l which is also given by Q and the result ^ is proved. □ 

3.4 Proof of Theorem 2, part 3 (local comparison around a = T) 

We fix V' ^ 1. We study the variations of ^^^^e^^" around a = T. From ([7]), we know: 
2Ko ' P(t) I \'^'p(t + a) 

therefore 

^Ap+Ko ^ ^ Pit) \ 

da 2Ko daY^ 'P{t + a)/' 
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Recalling that P depends on a (as N and Xp do), we have 

d dP 

-P[t + a) = -[t + a) + Pit + a). 

We then split the computations 

dXp + KQ.^ /,,,dP,, 1 \ Pit) (dP, , 

1 fdP,, P{t) dP , 
For a = T, the first term vanishes, and P{t + a) = P{t) i.e.. 

To compute this we again make use of the ODE (15|) which we multiply by — 
Averaging on a period we still get, for a = T, 

-^(*)\ _ \_ K-_/„;,2\ , nT^ /„i.2\^-\Fa 



Using dUl, we anive at 

We now have the derivative at a = T, 

^, ^^e^- = -AH(V'^) - 1). (9) 

We use here the notations A^(r) for ^| _^ and Xf(T) = Xp{T) to recall that we are 
studying the local behavior of Ai? and Xp around a = T,{'il) is fixed). We can directly compute 

da|„=T 2Ao 2Kq 2Kq 

Therefore, using (01) and Q, we obtain 



da\a=T 2Ko 



so that, using ([4]) and Q, we have 



Similarly we have 



Therefore, 



A'p(T) 



T + 



2Ko 
\p{T) 



T + 



\'p{T)-\'AT) 



2Ko 

MT){m - 1) 
r + 



2Kq 



Thanks to corollary |4l Ap (T) is positive. The assumption ([T]) leads to 

(^2)_l = ((^_l)2)>0. 

Finally we obtain 

A'p(r) - \'p{T) > 0, (10) 
and the second statement of the theorem follows then immediately from dJ) and (ITOl ). □ 

4 Modeling chronotherapy 

In the following we propose a model for chronotherapy by the introduction of a periodic death 
rate due to the effect of a drug on our cell division cycle model. 

4.1 Limit of single-phase division models 

We consider a population of cells following a general division equation with apoptosis rate d. 
As above, all coefficients are T-periodic with respect to time. 



f n(t, x) + £n{t, x) + {d{t, x) + K{t, x))n{t, x) 



0, 



i n{t,0) =2j^ K{t,x)n{t,x)dx. 
We consider the Floquet eigenproblem associated with this equation 

' §-^N{t, x) + ^Nit, x) + {d{t, x) + K{t, x) + XF)N{t, x) = 0, 
N{t, 0) = 2 K{t, x)N{t, x)dx, 

^ N>0, N{t, x)dxdt = 1 

We propose to model the effect of chronotherapy by adding a time T-periodic, age-independent 
death rate 7(t) representing the effect of a drug (for instance we may consider 7 proportional 
to the quantity of drug in the body). The cell population now follows the equation 

J f n(t, x) + ^n{t, x) + [d{t, x) + K{t, x) + -fit)]n{t, x) = 0, 

\ n(t,0) = 2 K{t,x)n{t,x)dx. 



The Floquet eigenproblem for this equation reads 

' |AfT(t, x) + |:iV^(t, x) + x) + K{t, x) + ^{t) + Xl)N''{t, x) = 0, 
< N"i{t, 0) = 2 /~ K{t, x)N"i{t, x)dx, 
^ >0, T - periodic j"^ N^{t, x)dxdt = 1. 
Lemma 5. The Floquet eigenvalue defined above satisfies 

Proof. We define 7 = 7— (7), T{t) = j{s)ds. Noticing that T is T-periodic, we define the 
function M by M{t, x) = N(t, x)e^^^\ It satisfies 

' |M(t, x) + -§^M{t, x) + {d{t, x) + K {t, x) + 7(t) + Af - (7)) 2;) = 0, 
< M{t, 0) = 2 (t, 2:)M(t, x)dx, 

M > 0, T-periodic. 

Therefore A^ = Ai? — (7) and up to a renormalization M = N'^. □ 
This result expresses that with such a simple model, chronotherapy is inefficient, since 
changing the moment of administration of a drug (in symbols, changing 7(t) into 'y(t + 9) 
where ^? is a real number) has no effect on the growth rate. In other words, in such one-phase 
models, this effects depends on (7). Only the total daily dose of the drug is relevant! 

4.2 Using multiphase models 

We now consider more realistic multiphase models. We use the additional ingredient that the 
real cell division cycle is multiphasic because of the existence of checkpoints between phases 
(mainly at the Gl/S and G2/M transitions) at which it can be an^ested if genome integrity is not 
preserved. We consider a cell cycle model with / phases where I > 1 (for instance / = 4 if we 
want to represent the classical phases G1-S-G2-M). We study / populations of cells, nj(t, x) 
being the density of cells of age x in phase i at time t. We use the convention 1 + 1 = 1 

^ni{t,x) + -§^ni{t,x) + [Ki_,i+i{t,x) + di{t,x)]ni{t,x) = 0, 
?ii+i(t,0) = Ki^i+i{t,y)ni{t,y)dy, 1 < i 

(1) 

ni(t,0) = 2 /q Kj^i{t,y)nj{t,y)dy, 
ni{0, x) = n^{x) given. 

Here represents the transition rate from phase itoi + l. At the end of phase / division 

occurs with rate Kj^i. To be as general as possible, we have considered death rates di in phase 
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i. As above, the coefficients are time T-periodic and we can consider the Floquet eigenproblem 



Ni{t,x) + i^Ni{t,x) + [Ki^i+i{t,x) + di{t,x) + X]Ni{t,x) = 0, 



Ni+i{t,0) = Ki^,+i{t,y)N,{t,y)dy, 1 < i 

Ni{t,0) = 2 Ki_,i{t,y)Nj{t,y)dy, 

Ni>0, T - periodic , li /o°° ^idxdt = 1. 



(2) 



We also consider the adjoint eigenproblem 



f d 



Mt,x) + -^(j)iit,x)-[Ki. 



n+1 



{t,x) + di{t,x) + X](j)i{t,x) 



§lh{t,x) + £(l)i{t,x) - [Ki_i{t,x)+di{t,x) + A]</./(t,x) = 2Ki^i(j)i{t,Q), 
(pi > 0, T — periodic, Yl fo° Nicpidxdt = 1. 

i 

(3) 

To model the effect of chronotherapy, we consider a cytotoxic drug acting only on a specific 
phase (for instance 5-Fluorouracil acts on S-phase, see IITtI for instance and the references 
therein) and, as in the previous section we represent its action by an additional death rate in 
phase j, 7(t) (we replace in phase j dj by dj + 7) . We also define eigenelements for the 
modified equation (A'^, N'^ , (fP). We multiply the first Une of (O (version with dj replaced by 
dj + 7, Ni by N] and A by A"^) by (pi, and © by N^. Summing over i and integrating over 
age and time, we obtain 



1 fOO 



, JO JO 



N/(j)idxdt 



7(t) / N](t)jdxdt 



(4) 



We shall not have here the problem encountered with one-phase models. We study the effect 
of a death rate ^{t + 6). We denote A^'^, N^'^ the eigenelements associated to an additional 
death rate £'^{t + 6) in phase j. We define 9) by 



F{e,9) = A - A- 



!le^{t + e)j^Nf(l>,dxdt 
Jo Jr Nf<t>idxdt 



(5) 



As we have A = A*^'^ for any 9, F(0, 9) = 0. Particularly it does not depend on 9. The question 
is: does F{e,9) depend on 9 for fixed e? To assess this question, we compute using dominated 
convergence 



dX{e,i 



\e=o= lim 



Fie, 9) 



j{t + 9) / Nj(l)jdxdt. 



(6) 



Therefore if neither the function 7 ( . ) nor the function Nj (j)j{.,x)dx are constant (contrarily 
to one-phase models, there are no compensating effect making Nj(j)j{., x)dx constant, see 

Fis 9) 

for instance the computations of the appendix), then lim — depends on 9 (we mean it 



e-»0 
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is not a constant function of 9) and so is (at least for small e) F{e, .). In this case the Taylor 
first order approximation around of A: A(e, 6) X + e ^{t + 9) Njcpjdxdt is not 
a constant function of 9 and neither is X{e,9), at least for small values of e. We illustrate 
this property numerically in the next section (see figure [5]l. It seems that the Taylor first order 
approximation is a very good approximation of the growth rate for a reasonable range of values 
of the amplitude e. 



5 Numerical simulations 

We illustrate the theorems proved above by several numerical simulations. We firstly present 
the numerical scheme, then we give several algorithmic properties. Finally tests are presented. 



5.1 Discretization 

In our numerical simulations we consider a pure division model : 

\ n{t,0) = 2Koi^{t) J^n{t,x)dx. 



0, 



(1) 



Consider time and age increments At, Ax and denote by k-i and ip^, the quantities ki = 
A'oX[a,+oo[(^^^) ^rid tp^ = 'iJj{kAt). Choosing first order finite differences, we obtain from 
equation ^ the following approximation with an error of order 0(| At| + | Ax|) 



n 



k+l 



+ 



0, 



Ki < L 



At Ax 

where {0 . . . /} is the set of all values of i to be considered in the discretization. Taking 
At = Ax (CFL = 1), we obtain the following compact discretization scheme: 



n. 



k+l 



1—1 



n, 



k+l 



2^" KjnfAt. 

0<i</ 



(2) 



Assume ip is periodic of period T > and consider a grid over [0, T] x [0, 1 At], consisting 
of squares with sides of length At = T/Nt, for some Nt G N (and / large enough, partic- 
ularly I At > a and / + 1 > Nt). Then, the populations at time {k + l)At for all ages in 
[0, 1 At] can be obtained from the corresponding populations at time A; At as follows: 



n 





fc+i 



V 



1 







2i/)'=K^_^ At 2tp''KjAt \ 





l+Ati/i^+ift^ 







h 



(3) 



It is clear that the matrix in ^ depends only on the time index k and is periodic of period 
A't^. We denote this matrix and the vectors respectively and n^'"*"^. The equation (|3]l 
can be written ri'^'+^ = Mj^n^. 
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5.2 Approximating the eigenvalue 

The algorithm has already been discussed in |[25l . We recall that the growth rate is defined as 
the unique real Xp such that ([T]) admits solutions of the form N{t,x)e^^* with > and 
N{.,x) is periodic. We can approximate it thanks to: 

Lemma 6 (Discrete Floquet theorem). 

There exists a unique real A and a unique sequence of vectors {M^)keN . ■^'^ 
such that 

1=0 

k 1-^ (A/"^) is NT-periodic , 
n'^, defined by = _A/^e^-''^* is solution to ^ 

Proof. The proof is standard and we recall it for the sake of completeness. It is based on the 
Perron Frobenius theorem. First we prove uniqueness. Supposing there exists such n'^, we 
have 

Mon°, 

Min^ = MiMon°, 

(7) 

Mku'' = MkMk-i . . . MiMon°, (8) 
MNT-iMNr-2 ■ ■ ■ MiMon°. (9) 

□ 

Lemma 7. The matrix M is nonnegative and primitive (and therefore is irreducible). 

Proof. The nonnegativity is obvious. To prove the primitivity, the key point is / + 1 > N-r and 
/At > o + 2At. For some e > we have for any k, if we denote by Idfc the identity matrix of 
order k, 

^ / 0. . .0 11 \ 

Notice that W is the Wielandt matrix of order / + 1 which is known to be primitive (see |[T5l ). 
Therefore for some p, W'^ > and thus for qNx > p, 

which yields the primitivity of M, the spectral radius of which, denoted here by p is then posi- 
tive. We denote by p its spectral radius. We have p > 0. 



0<i<I 



(4) 



(5) 
(6) 



We define 

thus, ^ reads n^^ = Mn°. 
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Back to the proof of the discrete Floquet theorem, we have 



Hence we have MnP = e^'^n^. This means that v? is a positive eigenvector of M associated 
to a positive eigenvalue e^"^ . From the Perron-Frobenius theorem, e'^'^ = p and nP = J\f^ is 
the (unique) associated eigenvector. The solution is unique. 

Conversely, if we know the Perron eigenvector V and the Perron eigenvalue p of M, then the 

/ ,\keN 
sequence ( A/" j defined by 

r A/'o = y, 

satisfies dH),® and ©for A = log(p(M)). □ 
For multiphase models, the idea is mainly the same. To compute p = e"^^ the spectral 
radius of M, the power algorithm is used. It converges thanks to the primitivity of M. 

5.3 Numerical results 

First we present some numerical results to illustrate theorem [2] We scale T = 1. We fix the 
value of Kq to 2 and test various periodic function ijj. We plot the curves 

a Air(a,V'), 
a — > Ap(a) 

We recall that the eigenvalues for the PeiTon problem can be directly computed thanks to 
lemma[3] From theorem|2l we know that these curves cross for x-coordinate a = T, the second 
part of the theorem teUs us that we expect (locally) the curve for Ai? to be above the curve for 
Ap for a < r and below it for a > T. We plot the curves A = Ap(a) and A = XF{a,Tp) for 
our functions if; and look at the crossing of curves around T (on the simulations, T = 1). We 
also give a more global view of \p{a,'il)sin) and Ap in figure |3] to illustrate the fact that the 
comparison is only local. Here, the parameters h and 8 are respectively set to 3 and 0.3. 



Name of the function 


Formulation on the interval [0, 1[ 




il)sq (square wave) 


l-8X[0,l/2[(t) + 0.1X[1/2,1[ 


1.81 


il}pk (peak function) 


0.1 + ht/5x[mi^) + {2h - ht/6)x[5,25[i^) 


1.99 


il)sin (sinusoidal) 


1 + 0.9cos(27rt) 


1.405 



Table 1 : Functions ^jj for the simulations 



From the last part of the demonstration of theorem [2j we expect, 

dXF{a,tppk) ^ dXF{a,ipsq) ^ dXF{a,tpsin) 
da da da 
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0.40 H : 1 : 1 : 1 : 1 : 1 : 1 : 1 : ' 

0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 

minimal age for division a 



Figure 2: Crossing of the Perron and Floquet curves (detail) for ijj = i/jsi 




0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 

minimal age for division a 



Figure 3: Crossing of the Perron and Floquet curves for ip = ipsi 



15 




Figure 4: Crossing of the Perron and Floquet curves for = tpsin (dash dot), ipsq (dots) and ippk 
(long dash). 

we give figure |4] as a confirmation. Finally we give some simulations to illustrate our remarks 
on chi-onotherapy. 

For the chronotherapy simulation we use the following parameter: we fix / = 3 (we 
consider S and G2 as a single phase). The parameter 7 is a periodic function (with strong 
variations on a period to have a stronger effect of the parameter 6 € (0, 1)). We compute the 
eigenvalue for a death rate in phase 2 (phase S-G2) having the value ejlt + 9). We test several 
value of e to determine whether or not the amplitude of the death rate changes the relative 
behavior of the eigenvalue with respect to O.The coefficients have the form: 

where Ki , a, are positive, ipi is a positive 1 periodic function. We give a simulation for the case 
described in the appendix (a case for which we can compute explicitly N24>2{t, x)dx). We 
fix Ki = 10 for all i, ai = 10/24, 02 = 12/24 = 0.5,03 = 2/24, ^(t) = 1 + 0.9cos(27rt) 
and ipi defined from ip as in the appendix. We choose ^{t) = cos^(27rt). With these choices 
of coefficients, we compute 



/>oo 

/ N2(t)2{t,x)dx = C -C'sm{2'Kt), 
Jo 



where C and C are positive constants. Therefore, 

\e,e _ \0 

lim = C + C' sin(27r0). 

In figure [51 we remark especially that the location of the optimal phase does not depend on 
e (since we have ^optimal = \ whatever the value of e) and con^esponds exactly to the value of 
9 maximizing siii(27r&), i.e., minimizing j{t + 9) N2(j)2dxdt. 
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Figure 5: Variation of the Floquet eigenvalue with respect to the parameter 6 for various amplitude 
for fixed 7 and amplitude e = 0.1, 0.5, 1 (from left to right). 

Concluding remarks 

The results of the present paper show that the periodic control on the transition rate Ki^i^i of 
cell cycle models yields richer behaviors than in the case in which only the death rates di are 
subject to a periodic control HI 13. In particular, the inequality of HI 13 does not caiTy over 
This is, to our knowledge, the first time that such results are shown -on special cases of the 
control- analytically, thus confirming numerical results first shown in |[8ll3- 

Our results also indicate that multiphase cell proliferation models are the simplest candi- 
dates to represent the effects of chronotherapy. Indeed, as shown in section |4~T1 in single-phase 
models, in the simple case when only death rates di are controlled by a periodic forcing term, 
the growth rate A is modified by a term depending only on the average over a period of the 
forcing term, so that no phase of the periodic control function can be relevant to account for 
differences in the resulting growth rate, contrary to what is observed in chronotherapy ifTSl . 
Furthermore such multiphase models take into account the existence of multiple checkpoints , 
and we know from cell cycle physiology that the minimal number of checkpoints to consider 
is 2: atGl/S and G2/M. 

We performed numerical and graphical results of section [51 on a 3-phase model with 1- 
periodic control on all phase transition functions Kj^j+i, where one represents chronotherapy 
as a 1 -periodic death term of amplitude e acting on the second phase (S/G2) only. These 
preliminary computational results, in particular performed in a simple analytically tractable 
case, seem to indicate that the effect of a chronotherapy on the growth rate X{e,9) highly 
depends on the amplitude e of the death rate but that the optimal phase 9 (related to the best 
peak infusion phase) is independent on e (see figure[5]l. In future work, we intend to introduce 
also an effect of chronotherapy on the ti'ansition rates Kj^j+i. 

From a more general point of view (i.e., independently of chronotherapeutic considera- 
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tions), Theorem |2] analytically shows, at least in the single-phase case, that under the control 
of a periodic function exerting its influence on cell division, a selective advantage is given 
to those cells that are able to divide with a cell cycle duration slightly lower than the control 
function period. But, as numerically illustrated on Fig.|3l this cell cycle duration should be not 
too much lower than the control period, or else the advantage is lost. This leads to a biological 
speculation (or prediction): in a population of proliferating cells with variable cycle duration 
times, all being under the control of a common 24 h-periodic circadian clock, those cells that 
are well controlled by the clock, and endowed with a cycle duration between say 2 1 h and 23 h 
should quickly outnumber the others. Hence in proliferating healthy tissues (fast renewing tis- 
sues such as gut or bone marrow), an intrinsic cell cycle time of 21 to 23 h should be observed 
(if such an observation is possible). 

Now to explain the initial tumour growth data that first motivated this study, we can spec- 
ulate in the following way: tumour cells are less sensitive than healthy cells to circadian clock 
control (indeed it is known from chronotherapeutics in oncology that "in contrast with consis- 
tent rhythmic changes in drug tolerability mechanisms in host tissues, tumour rhythms appear 
heterogeneous with regard to clock gene expression and rhythm in pharmacology determi- 
nants as a function of tumour type and stage" |[T8l ). so that their proUferation is more likely to 
be governed by a simple Perron eigenvalue rather than by one of the Floquet type. Tumour 
surrounding healthy cells and host immune cells, in contrast with tumour cells, ai^e still un- 
der circadian control and they may thus have a selective advantage over cancer cells as long 
as this circadian control is present. Circadian clock disruption by perturbed Ught-dark cycle 
destroys this advantage, and these perturbed host cells oppose in a less efficient way local tis- 
sue invasion by cancer cells, hence the resulting curves shown in the introduction. Of course 
such speculation remains to be documented (in particular by investigating differential circadian 
clock control on proliferation in tumour and healthy tissues), but this is our best explanation 
so far for this phenomenon. 

A Appendix 

A.l Existence theory for A 

This part is dedicated to the demonstration of the existence of the Floquet eigenvalue. Partic- 
ularly, we try to prove it under general hypothesis on the periodic function ijj. For instance, a 
short adaptation of the demonstration given in [i20il would be sufficient for the case of a pos- 
itive continuous periodic function ip, but one would like to have the possibility of studying 
non smooth functions such as a square wave (which for instance could have value 1 during the 
day and during the night). We give a proof of the existence of the Floquet eigenvalue in the 
one-phase model. It can easily be adapted for a multiphase-model with out death rates where 
the coefficients would have the form ETj^^j+i = Kiipi(t)x[a^^+oo[ with the same hypothesis 
on the functions ipi. We prove here existence of a solution to three eigenproblems: the direct 
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eigenproblem 

' dtN{t, x) + d^N{t, x) + (A + Ki;{t)x[a,+oo[{x))N{t, x) = 0, 

N{t, 0) = 2Ki:{t) N{t, x)dx, (10) 
^ TV > 0, Ndx = 1, 

the dual eigenproblem 

r -at</>(t, x) - + (A + i^V(t)X[a,+oo[(a;))</'(t, x) = 2Ki,{t)x[a,ooli^)^it^ 0)' 

1 > 0, Ncpdx = 1, 

(11) 

and the delay differential equation 

P{t) = -{Kil,{t) + X)P{t) + IKe-^^tpit - a)P{t - a), P > 0, [ P{t)dt = 1. (12) 

JO 

We give a normalization for P to ensure uniqueness. 

Theorem 8. For any positive T-periodic bounded function ^ 7^ 0, a > here exists a unique 
A, A^, (f), P such that P > is solution to ^ and N > is solution to ^ (N > if il^ is 
positive). 



The proof is based on the Kiein-Rutman theorem (see | fTO| | for instance). We consider a 
T-periodic nonnegative bounded function 0. We adapt the proof from |[20l to our case. 
First, using the methods of characteristics for the partial differential equations, we reduce the 
eigenproblems to integral equations on N{t,0), 0(i,O) and P. We consider three operators 
depending on a parameter ^u. For a bounded T-periodic function A4, we define A/i = Ci{A4) 
by 



Miit) = 2K iP{t-x)e-f'''-^f^^'-'-''+'^'^'Mit-x)dx, (13) 

J a 
poo 

N2{t) = 2K ij{t)e-''''-^ ^^i-^+'^d'Mit - x)dx, (14) 



/•oo 

M3{t) = 2K / il){t + x)e-^'''-^^>'^*+'^^'M{t + x)dx. (15) 
J a 

These operators are defined such that for /i = A, we get, 

p{t) = ciipm, 

iV(t,0) = £2(iV(.,0))(t), 

This means that the functions should be nonnegative eigenvectors of these three operators 
associated to the eigenvalue 1. 
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Lemma 9. For > 0, 

• Ci maps L^^.(0, T) into itself, 

• Ci,Ci are continuous compact operators on Cper(0, T), 

• Ci,Ci are strongly positive and £2 is nonnegative (strongly positive if ijj > Oj. 
Proof. For M bounded, one has, since for x > a, ■ip{t)dt < {ip){x — a) + \\il>\\ooT, 



\\Ci{M)\U < 2|^e™-^||Al|U = C\\M\U 

For continuity and compactness we only explicit the proof for i = 1, the case i = 3 is very 
similar. We consider A4 continuous and h small, 

POO 

Mi{t + h) = 2K iP{t + h-x)e-^'''-^^^^^'+^-''+'^'^'M{t + h-x)dx, 

J a 

/•oo ^ 
J a—h 

= 2K r iP{t-x)e-f'^''+^^-^^^^''^^^-''+'^'^'M{t-x)dx 

J a—h 

fOO 

J a 

= Ah + 2K ^{t - 2;) J^e-'^(^+'')-^^<^^'''^(*-^+')'^' - e-^^-^/<^''^(*-^+")'^"^7W(t - x)d 

POO 

+ 2K ilj{t-x)e-^'^-^-^>^*-^^'^'^'M{t-x)dx, 

J a 

= Ah + Bh+Mi{t). 
We have bounds on Ah and Bh, 

\Ah\ < 2K\\'4)\\\\M\\ooK 

\Bh\ < i^ll^lloo/illA/llloc < CKU\\oo\\M\\ooh. 

Therefore, using (lA-lKdATI ). we obtain the continuity and the compactness of operator Ci. 
Using the same techniques we can prove continuity and compactness of operator £3. The 
operator C2 needs regularity on ijj to be compact (and continuous). All these operators are 
positive. We can apply the Kiein-Rutman theorem (weak form lITOl ). We denote pi,P'i the 
spectral radii of respectively £1,^3. They are positive (since £(1) > e > 0, pi > e), so are 
the associated nonnegative eigenfunctions. If M.i{t) = 0, then 

il){t - x)Mi{t - x) = for x>a, 

which leads to ipM.i = and piM.i = O.Therefore Mi and similarly M.^ can not vanish. 

Lemma 10. 

Pi = P3- 



X 
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Proof. The first point is a straightforward computation. The second point uses the duality of 
operators L2 and £3, 

i:2[^M{){t)M?,m^ = r mMiit)C3{M3){t)dt, 

Jo 

Pi [ '^j{t)Mi{t)Mz{t)dt = P3 [ i){t)Mi{t)M3{t)dt. 
Jo Jo 

The existence of a solution to ^ is equivalent to the existence of a positive fixed point of £1 
for p = X, therefore, we need to find p such that pi{p) = 1. For ^ = 0, we have 

/:3(i) = 2. 

Therefore pi{0) = 2. As p is a decreasing function of p and /5i(oo) = 0, there exists some 
positive A such that pi{X) = 1. The solution to dSjl is then given by such A and P = Aii. 
Then, the function defined N{t, 0) = ip{t)A4i{t) and the characteristics 

N{t, x) = N{t - X, 0)e-^^-^" K^{t-x+s)xia.^i{s)ds ^ 
is solution to ©. We remark then, as A > 0, that N{t, 00) = 0. Similarly, we define (p by 

POO 

<t>{t, x) = / Ki,{t + y - x)x[aM^ym + y-x, 0)e- W(*+--)x[.,^[W'^^dy. 

J X 

This is a solution to ([TTI ). 



A.2 Explicit solutions for the multiphase eigenproblem 

In the following T = 1. 

We give here explicit solutions to the eigenproblem in the multiple phases case. We do not 
give details for the demonstration. We consider a 3 phase model without death terms, where 
the transition terms have the form: 

Ki^i+l{t,x) = KiiJi{t)x[ai,oo[{x)- 

Here, ipi is a positive 1-periodic function satisfying {ipi) = 1. We consider the following very 
specific case: we choose ai, 02, 03 > such that ai + 02 + 03 = 1, and we choose ipi in the 
following way, for a fixed positive 1-periodic function ip, 

Mt) = Ht), 

tp2{t) = il^{t-a2), 
fp3{t) = ip{t-a2-a3). 

To explain the form of the coefficients, we make the following remark: if we denote Pi{t) = 
Ni{t, x)dx (the same idea as for the one phase model), the 1-periodic functions Pi satisfies 
a system of delay differential equations and since oi + 02 + 03 = 1, the 1-periodic functions 
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Qi defined by Qi{t) = Pi{t), Q2{t) = P2{t + 02), Qi,{t) = Pz{t + 02 + 03) satisfy a system 
of ordinary differential equations. 

, I Qi{t) \ I -X-K,m IKse-^'^^m \ I Qi{t) \ 

- Q2it) ]=-[ Kie-^-^m -K2m-^ Q2{t) . 

V I \ K2e->^''^m -Ksm - A / V Qsii) J 

We denote M{t) the above matrix. Due to tlie special form of the functions tpi, we have 
M{t)M{t') = M{t')M{t), for all t, t'. Therefore we can write 

g(t) =exp(^^ M{s)ds^Q{0). 

The vector Q{t) is 1-periodic, thus, Q{0) has to be a positive eigenvector of exp( M{s)ds) 
associated to the eigenvalue 1. The matrix exp(Jj| M{s)ds) has eigenvalue 1 if and only if 

{Ki + X){K2 + X){K3 + A) - 2KiK2K3e~^^''^+''''+''-'^ = 0. (16) 

This leads to = e^^'('^W-^)'^''Qi(0), where Q(0) is apositive vector satisfying M{s)dsQ{0) = 
0. Then, we can compute Pi{t) and Ni{t, 0). Finally, using the methods of characteristics, the 
eigenfunctions Ni are given, up to a normalization, by 

N2{t,x) = KiUiip{t - x)e^'^o"''(V'(«)-i)'^«-^^-/o"^2»/'(i-2:+s-a2)x[a2,oo[(s)rfi*^ 
iV3(t,x) = K2U2ipit - 02 - x)e^^o~""'''('^W-i)'^^-^^-/o'^3^(*-^+^-«2-«3)x[a3,^[(«)rf^^ 
where 







( 


1 \ 


U2 






ii'2+A 






[ 


2X36-^"! / 



The adjoint eigenfunctions are given by the formulas 



where 

M^2 = K^e-^^i 
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Basically, the ideas for the computations of (pi are the same, based on the following remark, as 

/•CO 

Ut,x) = / i^,^i+i(t + y,x + y)<A,,4-i(i + y,0)e-/ci'A+i^-»+i(t+y',-H-yWrfy, 



(with a factor 2 for i = 3), we have (j)i{t,x) = <j)(t,ai) for a > ai. This leads to a 
differential equation for (f)i{t,0). Details are left to the reader. In this case, we compute 
/o°° ^iit, x)(j)i{t, x)dx. As we have x) = (j)i{t, ai) for x > ai, 

f-ii r-oo 

Ni(t,x)(f)i{t,x)dx = / Ni{t,x)cl)i(t,x)dx + (j)i{t,ai) / Ni{t,x)dx. 







We have 



Ni{t,x)Mt,x) = iKi + X)e^''' / Tp{t - X + ai)dx + UiVie^''^ 







N2{t,x)4>2{t,x)dx = (Ki + A)e^"i / i;{t - x)dx + U2V2e^''\ 



N3{t,x)(l)3it,x)dx = (Ki + A)e^"i / 7p{t - a2 - x)dx + UsVse^''^ 

Jo 

Particularly, in this case, Ni(j)dx is not always constant. We denote ^{t) = jQ{ilj{s) — l)ds, 
it is a 1 periodic function. We also denote Cj = UiVie^'^\C = {Ki + A)e^"\ both these 
constants are positive, 

/•oo 

/ Ni{t,x)(l)i{t,x) = C7(ai + ^(t)-^'(t + ai)) + C7i, 
Jo 

/■oo 

/ N2{t,x)(l)2{t,x)dx = C7(a2 + ^(t-a2) -^(t)) + C72, 

N3{t,x)(l)3{t,x)dx = C(a3 + ^'(t + ai)-^'(t-a2)) + C3. 

For instance, using the parameters of the simulation, we have, ^(t) = ^ sin(27rt), 02 = 0.5, 
f°° 9 

/ N2(t)2{t, x)dx = {Ca2 + C2) - 2C— sin(27rt) = C" - sin(27rt), 
Jo 27r 

f'\ POO /•! /•! 

/ jit+e) N2(t>2{t,x)dx = C' I cos^{2T:{t+e))dt-C2 I cos^(27r(t+6')) sin(27rt)(it, 
Jo Jo Jo Jo 

a short computation leads to 

7(t) = — cosfGvrt) H cos(47rt) H cos(27ri) H , 

^ 32 ^ ^ 24 ^ ^ 32 ^ ^ 16 

therefore, 

/•I /•oc 

/ 7(t + e) / N2(p2{t,x)dx = C" - C2 sm{27re), 
Jo Jo 
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Where C" > > 0. Therefore, in this particular case, 

\e,d _ \0 

lim ^ = C'{ sin(27r0) - C" . 
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